
#### Maps Venezuela ### 

library(tidyverse)
library(readxl)
library(stringi)
library(sf)
library(ggtext)
library(grid)
library(ggpubr)

vene <- st_read("data/Venezuela_municipio_capas_vectoriales/ba_4189_dpt_municipal_de_venezuela_250000_20150707.shp")

# Convert specific columns to UTF-8
vene <- vene %>%
  mutate(
    municipio = iconv(municipio, from = "UTF-8", to = "ISO-8859-1"),
    cap_mun = iconv(cap_mun, from = "UTF-8", to = "ISO-8859-1"),
    ent_fdral = iconv(ent_fdral, from = "UTF-8", to = "ISO-8859-1")
  )

# Define custom replacements for accented characters
remove_accents <- function(text) {
  text %>%
    str_replace_all(c(
      "á" = "a", "é" = "e", "í" = "i", "ó" = "o", "ú" = "u",
      "Á" = "A", "É" = "E", "Í" = "I", "Ó" = "O", "Ú" = "U"
    ))
}

vene <- vene %>% mutate(
  cod_ent = case_when(
    ent_fdral == "DISTRITO CAPITAL" ~ "01",
    ent_fdral == "ESTADO AMAZONAS" ~ "22",
    ent_fdral == "ESTADO ANZOÁTEGUI" ~ "02",
    ent_fdral == "ESTADO APURE" ~ "03",
    ent_fdral == "ESTADO ARAGUA" ~ "04",
    ent_fdral == "ESTADO BARINAS" ~ "05",
    ent_fdral == "ESTADO BOLÍVAR" ~ "06",
    ent_fdral == "ESTADO CARABOBO" ~ "07",
    ent_fdral == "ESTADO COJEDES" ~ "08",
    ent_fdral == "ESTADO DELTA AMACURO" ~ "23",
    ent_fdral == "ESTADO FALCÓN" ~ "09",
    ent_fdral == "ESTADO GUÁRICO" ~ "10",
    ent_fdral == "ESTADO LARA" ~ "11",
    ent_fdral == "ESTADO MÉRIDA" ~ "12",
    ent_fdral == "ESTADO BOLIVARIANO MIRANDA" ~ "13",
    ent_fdral == "ESTADO MONAGAS" ~ "14",
    ent_fdral == "ESTADO NUEVA ESPARTA" ~ "15",
    ent_fdral == "ESTADO PORTUGUESA" ~ "16",
    ent_fdral == "ESTADO SUCRE" ~ "17",
    ent_fdral == "ESTADO TÁCHIRA" ~ "18",
    ent_fdral == "ESTADO TRUJILLO" ~ "19",
    ent_fdral == "ESTADO YARACUY" ~ "20",
    ent_fdral == "ESTADO ZULIA" ~ "21",
    ent_fdral == "ESTADO VARGAS" ~ "24",
    TRUE ~ cod_ent  # Keep the original if no match is found
  ),
  
  # Ensure all codes are two digits
  cod_ent = str_pad(cod_ent, width = 2, pad = "0"),
  cod_ent_mun = paste(cod_ent, cod_mun, sep = ""),
  municipio = remove_accents(municipio),
  municipio = tolower(municipio),
  municipio = str_remove_all(municipio, "municipio turistico |municipio bolivariano |municipio |indigena bolivariano ")
)

mun_names <- vene %>% as_tibble() %>% select(municipio) %>% filter(!str_detect(municipio, "^zona ")) %>% distinct()
mun_names <- mun_names$municipio


ven_mun <- readRDS("ven_mun_elec_2006_2024_final.rds")

ven_mun <- ven_mun %>%
  mutate(
    cod_edo = str_pad(cod_edo, width = 2, pad = "0"),
    cod_mun = str_pad(cod_mun, width = 2, pad = "0")
  )

mun_elec_names <- unique(ven_mun$municipio)

mun_diff <- setdiff(mun_names, mun_elec_names)

mun_diff

vene <- vene %>%
  mutate(
    ent_fdral = case_when(
      ent_fdral == "DISTRITO CAPITAL" ~ "DISTRITO CAPITAL",
      ent_fdral == "ESTADO AMAZONAS" ~ "AMAZONAS",
      ent_fdral == "ESTADO ANZOÁTEGUI" ~ "ANZOATEGUI",
      ent_fdral == "ESTADO APURE" ~ "APURE",
      ent_fdral == "ESTADO ARAGUA" ~ "ARAGUA",
      ent_fdral == "ESTADO BARINAS" ~ "BARINAS",
      ent_fdral == "ESTADO BOLÍVAR" ~ "BOLIVAR",
      ent_fdral == "ESTADO CARABOBO" ~ "CARABOBO",
      ent_fdral == "ESTADO COJEDES" ~ "COJEDES",
      ent_fdral == "ESTADO DELTA AMACURO" ~ "DELTA AMACURO",
      ent_fdral == "ESTADO FALCÓN" ~ "FALCON",
      ent_fdral == "ESTADO GUÁRICO" ~ "GUARICO",
      ent_fdral == "ESTADO LARA" ~ "LARA",
      ent_fdral == "ESTADO MÉRIDA" ~ "MERIDA",
      ent_fdral == "ESTADO BOLIVARIANO MIRANDA" ~ "MIRANDA",
      ent_fdral == "ESTADO MONAGAS" ~ "MONAGAS",
      ent_fdral == "ESTADO NUEVA ESPARTA" ~ "NUEVA ESPARTA",
      ent_fdral == "ESTADO PORTUGUESA" ~ "PORTUGUESA",
      ent_fdral == "ESTADO SUCRE" ~ "SUCRE",
      ent_fdral == "ESTADO TÁCHIRA" ~ "TACHIRA",
      ent_fdral == "ESTADO TRUJILLO" ~ "TRUJILLO",
      ent_fdral == "ESTADO YARACUY" ~ "YARACUY",
      ent_fdral == "ESTADO ZULIA" ~ "ZULIA",
      ent_fdral == "ESTADO VARGAS" ~ "LA GUAIRA",
      TRUE ~ ent_fdral  # Keep original if not matched
    ),
    municipio = case_when(
      municipio == "las mercedes" ~ "juan jose rondon",
      municipio == "pedro zaraza" ~ "zaraza",
      municipio == "moron" ~ "moran",
      municipio == "francisco del carmen carvajal" ~ "carvajal",
      municipio == "pedro maria freites" ~ "freites",
      municipio == "miguel peña" ~ "valencia",
      municipio == "jose felipe maquez cañizales" ~ "jose felipe marquez cañizales",
      municipio == "san jose de guanipa" ~ "guanipa",
      municipio == "sir arthur mc gregor" ~ "mc gregor",
      municipio == "tinaquillo" ~ "autonomo tinaquillo",
      municipio == "tomas lander" ~ "lander",
      municipio == "fernando de peñalver" ~ "peñalver",
      municipio == "jose antonio paez" & ent_fdral == "BARINAS" ~ "pedraza",
      municipio == "simon bolivar" & ent_fdral == "BOLIVAR" ~ "bolivar",
      municipio == "simon bolivar" & ent_fdral == "ANZOATEGUI" ~ "bolivar",
      municipio == "francisco de miranda" & ent_fdral == "ANZOATEGUI" ~ "miranda",
      municipio == "ezequiel zamora" & ent_fdral == "COJEDES" ~ "san carlos",
      municipio == "francisco de miranda" & ent_fdral == "GUARICO" ~ "miranda",
      municipio == "jose felix ribas" & ent_fdral == "GUARICO" ~ "ribas",
      TRUE ~ municipio  # Keep all other values unchanged
    ),
    mun_edo = paste(ent_fdral, municipio, sep = "_")
      )

#### Merge wider dataset #####

ven_mun_wide <- ven_mun %>%
  select(year, estado, municipio, cod_edo, cod_mun, of_p, op_p,otro_p, turnout) %>% 
  pivot_wider(
    names_from = year,
    values_from = c(of_p, op_p, otro_p, turnout),
    names_glue = "{.value}_{year}"
  ) %>% 
  unite(mun_edo, estado, municipio, sep = "_", remove = T)

ven_elec_mapa <- vene %>% 
  left_join(., ven_mun_wide)

#write.csv(ven_elec_mapa, "ven_elec_mun_cart.csv", row.names = F)
#saveRDS(ven_elec_mapa, "ven_elec_mun_cart.rds")

#######################################

## Map plots ####

ven_elec_mapa <- readRDS("ven_elec_mun_cart.rds")

## Create a separate layer for Esequibo Territory (to overlay hatching lines) ####
esequibo_layer <- ven_elec_mapa %>% filter(entidad == "TERRITORIO ESEQUIBO") %>% 
  mutate(esequibo = "Reclamation")

names(ven_elec_mapa)


## Strongholds

library(RColorBrewer)

ven_elec_mapa <- ven_elec_mapa %>%
  rowwise() %>%
  mutate(
    # Compute mean and median vote share for each party across elections
    mean_of_p = mean(c(of_p_2006, of_p_2012, of_p_2013, of_p_2018, of_p_2024), na.rm = TRUE),
    mean_op_p = mean(c(op_p_2006, op_p_2012, op_p_2013, op_p_2024), na.rm = TRUE),
    mean_otro_p = mean(c(otro_p_2006, otro_p_2012, otro_p_2013, otro_p_2018, otro_p_2024), na.rm = TRUE),
    mean_turnout = mean(c(turnout_2006, turnout_2012, turnout_2013, turnout_2018, turnout_2024), na.rm = T),
    
    median_of_p = median(c(of_p_2006, of_p_2012, of_p_2013, of_p_2018, of_p_2024), na.rm = TRUE),
    median_op_p = median(c(op_p_2006, op_p_2012, op_p_2013, op_p_2024), na.rm = TRUE),
    median_otro_p = median(c(otro_p_2006, otro_p_2012, otro_p_2013, otro_p_2018, otro_p_2024), na.rm = TRUE),
    
    # Stronghold classification based on historical dominance (mean vote share)
    stronghold = case_when(
      mean_of_p > 55 ~ "Chavismo Stronghold",
      mean_op_p > 55 ~ "Opposition Stronghold",
      mean_otro_p > 55 ~ "Other Stronghold",
      TRUE ~ "Battleground"
    )
  ) %>%
  ungroup()

strongholds <- ggplot(ven_elec_mapa) +
  geom_sf(aes(fill = stronghold), color = "black", lwd = 0.1) +
  # Overlay dashed lines specifically over Esequibo territory
  geom_sf(data = esequibo_layer, aes(fill = esequibo), linetype = "dashed", size = 1) +
  scale_fill_manual(
    name = "Political Stronghold",
    breaks = c("Chavismo Stronghold",
      "Opposition Stronghold",
      "Other Stronghold",
      "Battleground"),
    values = c(
      "Chavismo Stronghold" = "red",
      "Opposition Stronghold" = "blue",
      "Other Stronghold" = "purple",
      "Battleground" = "gold"
    )
  ) +
  theme(line = element_blank(),
        axis.text = element_blank(),
        axis.title = element_blank(),
        panel.background = element_blank(),
        plot.caption = element_text(hjust = 0, face = "italic"),
                legend.text = element_text(size = 10),
                title = element_text(size = 14)) +
  labs(
    title = "Figure 2: Strongholds by Municipality",
    subtitle = "Average vote share: Venezuelan presidential elections 2006-2024",
    caption = "Note: Strongholds are municipalities where a party has historically received >55% of the vote.\nMain opposition vote share in 2018 is excluded from the average calculation due to non-participation.\nWhite municipalities correspond to missing data. Grey correspond to the Esequibo territory \nwhich is under reclamation."
  )
  
strongholds

ggsave("hist_strongholds.jpg", strongholds, width = 10, height = 8)

#### Average turnout #####

ven_elec_mapa <- ven_elec_mapa %>% 
  mutate(mean_turnout = case_when(
    is.na(mean_turnout) ~ lag(mean_turnout),
    TRUE ~ mean_turnout
  ))

map_turnout <- ggplot(ven_elec_mapa) + 
  geom_sf(aes(fill = mean_turnout)) +
  geom_sf(data = ven_elec_mapa, fill = NA, color = "black", lwd = 0.15) +
  geom_sf(data = esequibo_layer, fill = "darkgrey", linetype = "dashed", size = 1) +
  theme(line = element_blank(),
        axis.text = element_blank(),
        axis.title = element_blank(),
        panel.background = element_blank(),
        legend.text = element_text(size = 10),
        title = element_text(size = 14)) +
  guides(fill = guide_colorbar(title = "Participation rate (%)")) +  # Colorbar guide for continuous scale
  scale_fill_gradient(
    name = "Turnout \n",
    na.value = "white",
    low = "yellow", 
    high = "darkgreen",
    limits = c(40, 80), # Ensures correct range
    breaks = seq(40, 100, by = 10),  # Adjusts legend steps from 55% to 100%
    labels = function(x) paste0(x, "%")  # Format labels as percentages
  ) +
  labs(
    title = "Figure 1: Average electoral turnout by Municipality",
    subtitle = "Venezuela Presidential elections 2006-2024"
  )

map_turnout

summary(ven_elec_mapa$mean_turnout)

ggsave("map_turnout_2006_2024.jpg", map_turnout, width = 10, height = 8)

ven_state_turnout <- ven_elec_mapa %>%
  st_drop_geometry() %>%
  group_by(ent_fdral) %>%
  summarise(mean_turnout = mean(mean_turnout, na.rm = TRUE))

ven_state_map <- st_read("data/Venezuela_capasvectoriales_estado/ba_4189_dpt_estadal_de_venezuela_250000_20150707.shp")
  
  
  
  select(ent_fdral, geometry) %>%
  left_join(ven_state_turnout, by = "ent_fdral")



### Turnout 2013 ####

ven_elec_mapa <- ven_elec_mapa %>% 
  mutate(turnout_2013 = case_when(
    is.na(turnout_2013) ~ lag(turnout_2013),
    TRUE ~ turnout_2013
  ))

map_turnout_13 <- ggplot(ven_elec_mapa) + 
  geom_sf(aes(fill = turnout_2013)) +
  geom_sf(data = ven_elec_mapa, fill = NA, color = "black", lwd = 0.15) +
  geom_sf(data = esequibo_layer, fill = "darkgrey", linetype = "dashed", size = 1) +
  theme(line = element_blank(),
        axis.text = element_blank(),
        axis.title = element_blank(),
        panel.background = element_blank(),
        title = element_text(size = 12)) +
  guides(fill = guide_colorbar(title = "Participation rate (%)")) +  # Colorbar guide for continuous scale
  scale_fill_gradient(
    name = "Turnout \n",
    na.value = "white",
    low = "yellow", 
    high = "darkgreen",
    limits = c(50, 90), # Ensures correct range
    breaks = seq(50, 100, by = 10),  # Adjusts legend steps from 55% to 100%
    labels = function(x) paste0(x, "%")  # Format labels as percentages
  ) +
  labs(
    title = "Electoral Turnout in Venezuela",
    subtitle = "2013 Presidential Election"
  )

map_turnout_13

ggsave("map_turnout_2013.jpg", map_turnout_13)

library(cowplot)

ggarrange(map_turnout_13, strongholds, 
          nrow = 2) %>%
  annotate_figure(., top = text_grob("Figure 1. Electoral turnout and political strongholds in Venezuela 2006-2024",
                  face =  "bold", size = 14)) %>%  
  ggexport(filename = "maps_turnout_stronholds.jpg", width = 550, height = 750)

map_chavismo_13 <- ggplot(ven_elec_mapa) + 
  geom_sf(aes(fill = of_p_2013*100)) +
  geom_sf(data = ven_elec_mapa, fill = NA, color = "black", lwd = 0.15) +
  theme(line = element_blank(),
        axis.text = element_blank(),
        axis.title = element_blank(),
        panel.background = element_blank()) +
  guides(fill = guide_legend(title = "Chavismo Vote (%)")) +
  scale_fill_gradient2("Chavismo Vote \n", low = "white", mid = "yellow", high = "red", 
                       midpoint = 50) +
  labs(title = "Chavismo Vote Share in Venezuela",
       subtitle = "2013 Presidential Election")

map_chavismo_13

ggsave("map_chavismo_2013.jpg", map_chavismo_13)

map_opposition_13 <- ggplot(ven_elec_mapa) + 
  geom_sf(aes(fill = op_p_2013*100)) +
  geom_sf(data = ven_elec_mapa, fill = NA, color = "black", lwd = 0.15) +
  theme(line = element_blank(),
        axis.text = element_blank(),
        axis.title = element_blank(),
        panel.background = element_blank()) +
  guides(fill = guide_legend(title = "Opposition Vote (%)")) +
  scale_fill_gradientn(colours=(brewer.pal(n = 5, "Blues"))) +
  labs(title = "Opposition Vote Share in Venezuela",
       subtitle = "2013 Presidential Election")

map_opposition_13

ggsave("map_opposition_2013.svg", map_opposition_13)

ven_elec_mapa <- ven_elec_mapa %>%
  mutate(winner_13 = ifelse(of_p_2013 > op_p_2013, "Chavismo", "Opposition"))

table(ven_elec_mapa$entidad)


# Base map
map_winner_13 <- ggplot() + 
  # Electoral winners by municipality
  geom_sf(data = ven_elec_mapa, aes(fill = winner_13), color = "black", lwd = 0.15) +
  
  # Overlay dashed lines specifically over Esequibo territory
  geom_sf(data = esequibo_layer, aes(fill = esequibo), linetype = "dashed", size = 1) +
  
  # Map aesthetics
  theme_minimal() +
  theme(
    line = element_blank(),
    axis.text = element_blank(),
    axis.title = element_blank(),
    panel.background = element_blank()
  ) +
  
  # Custom fill colors
  scale_fill_manual("Winner \n", 
                    breaks = c("Chavismo", "Opposition"),
                    values = c("Chavismo" = "red", "Opposition" = "blue", "Reclamation" = "grey"),
                    na.translate = FALSE) +
  
  # Title and subtitle
  labs(
    title = "Electoral Winner by Municipality",
    subtitle = "2013 Presidential Election",
    caption = "Note: White municipalities correspond to missing data. Grey correspond to the Esequibo territory \nwhich is under reclamation."
  ) +
  
  # Formatting caption
  theme(plot.caption = element_text(hjust = 0, face = "italic"))

map_winner_13

ggsave("map_winner_2013.jpg", map_winner_13)


ven_elec_mapa <- ven_elec_mapa %>%
  mutate(winner_24 = ifelse(of_p_2024 > op_p_2024, "Chavismo", "Opposition"))

map_winner_24 <- ggplot(ven_elec_mapa) + 
  geom_sf(aes(fill = winner_24)) +
  geom_sf(data = ven_elec_mapa, fill = NA, color = "black", lwd = 0.15) +
  theme(line = element_blank(),
        axis.text = element_blank(),
        axis.title = element_blank(),
        panel.background = element_blank()) +
  scale_fill_manual("Winner \n", values = c("Chavismo" = "red", "Opposition" = "blue")) +
  labs(title = "Electoral Winner by Municipality",
       subtitle = "2024 Presidential Election")

map_winner_24

ggsave("map_winner_2024.svg", map_winner_24)


ven_elec_mapa <- ven_elec_mapa %>%
  mutate(
    of_dif_18_24 = (of_p_2024 - of_p_2018) * 100,  # Chavismo difference
    op_dif_18_24 = (op_p_2024 - op_p_2018) * 100,   # Opposition difference
    of_dif_18_24_cat = case_when(
      of_dif_18_24 <= -75 ~ "> 75% Loss",
      of_dif_18_24 <= -50 ~ "> 50% Loss",
      of_dif_18_24 <= -25 ~ "> 25% Loss",
      of_dif_18_24 <= -10 ~ "> 10% Loss",
      of_dif_18_24 >= 10 ~ "> 10% Gain",
      of_dif_18_24 >= 25 ~ "> 25% Gain",
      of_dif_18_24 >= 50 ~ "> 50% Gain",
      of_dif_18_24 >= 75 ~ "> 75% Gain",
      is.na(of_dif_18_24) ~ "NA",
      TRUE ~ "No Significant Change"
    ))

map_of_gains_24 <- ggplot(ven_elec_mapa) + 
  geom_sf(aes(fill = of_dif_18_24_cat)) +
  geom_sf(data = ven_elec_mapa, fill = NA, color = "black", lwd = 0.15) +
  theme(line = element_blank(),
        axis.text = element_blank(),
        axis.title = element_blank(),
        panel.background = element_blank()) +
  scale_fill_manual(
    "Chavismo Vote Change \n(2018-2024)", 
    values = c(
      "> 75% Loss" = "darkred",
      "> 50% Loss" = "firebrick2",
      "> 25% Loss" = "red",
      "> 10% Loss" = "black",
      "NA" = "grey",
      "No Significant Change" = "white",
      "> 10% Gain" = "pink",
      "> 25% Gain" = "red",
      "> 50% Gain" = "darkred",
      "> 75% Gain" = "black"
    )
  ) +
  labs(
    title = "Chavismo Vote Gains/Losses by Municipality",
    subtitle = "2018-2024 Presidential Election"
  )

map_of_gains_24

ggsave("map_of_loses_2024.svg", map_of_gains_24)




